Prony's method
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
top
Prony analysis (Prony's method) was developed by Gaspard Riche de Prony in 1795. However, practical use of the method awaited the digital computer.cite-ref-1[1] Similar to the Fourier transform, Prony's method extracts valuable information from a uniformly sampled signal and builds a series of damped complex exponentials or damped sinusoids. This allows the estimation of frequency, amplitude, phase and damping components of a signal.
Contents
โข The method
โข Representations
โข See also
โข Notes
โข References
โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ
The method
Let f ( t ) {\displaystyle f(t)} be a signal consisting of N {\displaystyle N} evenly spaced samples. Prony's method fits a function
f ^ ( t ) = โ i = 1 N A i e ฯ i t cos โก ( ฯ i t + ฯ i ) {\displaystyle {\hat {f}}(t)=\sum _{i=1}^{N}A_{i}e^{\sigma _{i}t}\cos(\omega _{i}t+\phi _{i})}
to the observed f ( t ) {\displaystyle f(t)} . After some manipulation utilizing Euler's formula, the following result is obtained, which allows more direct computation of terms:
f ^ ( t ) = โ i = 1 N 1 2 A i ( e j ฯ i e ฮป i + t + e โ j ฯ i e ฮป i โ t ) , {\displaystyle {\begin{aligned}{\hat {f}}(t)=\sum _{i=1}^{N}{\tfrac {1}{2}}A_{i}\left(e^{j\phi _{i}}e^{\lambda _{i}^{+}t}+e^{-j\phi _{i}}e^{\lambda _{i}^{-}t}\right),\end{aligned}}}
where
ฮป i ยฑ = ฯ i ยฑ j ฯ i {\displaystyle \lambda _{i}^{\pm }=\sigma _{i}\pm j\omega _{i}} are the eigenvalues of the system,
ฯ i = โ ฯ 0 , i ฮพ i {\displaystyle \sigma _{i}=-\omega _{0,i}\xi _{i}} are the damping components,
ฯ i = ฯ 0 , i 1 โ ฮพ i 2 {\displaystyle \omega _{i}=\omega _{0,i}{\sqrt {1-\xi _{i}^{2}}}} are the angular-frequency components,
ฯ i {\displaystyle \phi _{i}} are the phase components,
A i {\displaystyle A_{i}} are the amplitude components of the series,
j {\displaystyle j} is the imaginary unit ( j 2 = โ 1 {\displaystyle j^{2}=-1} ).
Representations
Prony's method is essentially a decomposition of a signal with M {\displaystyle M} complex exponentials via the following process:
Regularly sample f ^ ( t ) {\displaystyle {\hat {f}}(t)} so that the n {\displaystyle n} -th of N {\displaystyle N} samples may be written as
F n = f ^ ( ฮ t n ) = โ m = 1 M B m e ฮป m ฮ t n , n = 0 , โฆ , N โ 1. {\displaystyle F_{n}={\hat {f}}(\Delta _{t}n)=\sum _{m=1}^{M}\mathrm {B} _{m}e^{\lambda _{m}\Delta _{t}n},\quad n=0,\dots ,N-1.}
If f ^ ( t ) {\displaystyle {\hat {f}}(t)} happens to consist of damped sinusoids, then there will be pairs of complex exponentials such that
B i ( 1 ) = 1 2 A i e ฯ i j , B i ( 2 ) = 1 2 A i e โ ฯ i j , ฮป i ( 1 ) = ฯ i + j ฯ i , ฮป i ( 2 ) = ฯ i โ j ฯ i , {\displaystyle {\begin{aligned}\mathrm {B} _{i}^{(1)}&={\tfrac {1}{2}}A_{i}e^{\phi _{i}j},\\\mathrm {B} _{i}^{(2)}&={\tfrac {1}{2}}A_{i}e^{-\phi _{i}j},\\\lambda _{i}^{(1)}&=\sigma _{i}+j\omega _{i},\\\lambda _{i}^{(2)}&=\sigma _{i}-j\omega _{i},\end{aligned}}}
where
B i ( 1 ) e ฮป i ( 1 ) t + B i ( 2 ) e ฮป i ( 2 ) t = 1 2 A i e ฯ i j e ( ฯ i + j ฯ i ) t + 1 2 A i e โ ฯ i j e ( ฯ i โ j ฯ i ) t = A i e ฯ i t cos โก ( ฯ i t + ฯ i ) . {\displaystyle {\begin{aligned}\mathrm {B} _{i}^{(1)}e^{\lambda _{i}^{(1)}t}+\mathrm {B} _{i}^{(2)}e^{\lambda _{i}^{(2)}t}&={\tfrac {1}{2}}A_{i}e^{\phi _{i}j}e^{(\sigma _{i}+j\omega _{i})t}+{\tfrac {1}{2}}A_{i}e^{-\phi _{i}j}e^{(\sigma _{i}-j\omega _{i})t}\\&=A_{i}e^{\sigma _{i}t}\cos(\omega _{i}t+\phi _{i}).\end{aligned}}}
Because the summation of complex exponentials is the homogeneous solution to a linear difference equation, the following difference equation will exist:
f ^ ( ฮ t n ) = โ m = 1 M f ^ [ ฮ t ( n โ m ) ] P m , n = M , โฆ , N โ 1. {\displaystyle {\hat {f}}(\Delta _{t}n)=\sum _{m=1}^{M}{\hat {f}}[\Delta _{t}(n-m)]P_{m},\quad n=M,\dots ,N-1.}
The key to Prony's Method is that the coefficients in the difference equation are related to the following polynomial:
z M โ P 1 z M โ 1 โ โฏ โ P M = โ m = 1 M ( z โ e ฮป m ) . {\displaystyle z^{M}-P_{1}z^{M-1}-\dots -P_{M}=\prod _{m=1}^{M}\left(z-e^{\lambda _{m}}\right).}
These facts lead to the following three steps within Prony's method:
1) Construct and solve the matrix equation for the P m {\displaystyle P_{m}} values:
[ F M โฎ F N โ 1 ] = [ F M โ 1 โฆ F 0 โฎ โฑ โฎ F N โ 2 โฆ F N โ M โ 1 ] [ P 1 โฎ P M ] . {\displaystyle {\begin{bmatrix}F_{M}\\\vdots \\F_{N-1}\end{bmatrix}}={\begin{bmatrix}F_{M-1}&\dots &F_{0}\\\vdots &\ddots &\vdots \\F_{N-2}&\dots &F_{N-M-1}\end{bmatrix}}{\begin{bmatrix}P_{1}\\\vdots \\P_{M}\end{bmatrix}}.}
Note that if N โ 2 M {\displaystyle N\neq 2M} , a generalized matrix inverse may be needed to find the values P m {\displaystyle P_{m}} .
2) After finding the P m {\displaystyle P_{m}} values, find the roots (numerically if necessary) of the polynomial
z M โ P 1 z M โ 1 โ โฏ โ P M . {\displaystyle z^{M}-P_{1}z^{M-1}-\dots -P_{M}.}
The m {\displaystyle m} -th root of this polynomial will be equal to e ฮป m {\displaystyle e^{\lambda _{m}}} .
3) With the e ฮป m {\displaystyle e^{\lambda _{m}}} values, the F n {\displaystyle F_{n}} values are part of a system of linear equations that may be used to solve for the B m {\displaystyle \mathrm {B} _{m}} values:
[ F k 1 โฎ F k M ] = [ ( e ฮป 1 ) k 1 โฆ ( e ฮป M ) k 1 โฎ โฑ โฎ ( e ฮป 1 ) k M โฆ ( e ฮป M ) k M ] [ B 1 โฎ B M ] , {\displaystyle {\begin{bmatrix}F_{k_{1}}\\\vdots \\F_{k_{M}}\end{bmatrix}}={\begin{bmatrix}(e^{\lambda _{1}})^{k_{1}}&\dots &(e^{\lambda _{M}})^{k_{1}}\\\vdots &\ddots &\vdots \\(e^{\lambda _{1}})^{k_{M}}&\dots &(e^{\lambda _{M}})^{k_{M}}\end{bmatrix}}{\begin{bmatrix}\mathrm {B} _{1}\\\vdots \\\mathrm {B} _{M}\end{bmatrix}},}
where M {\displaystyle M} unique values k i {\displaystyle k_{i}} are used. It is possible to use a generalized matrix inverse if more than M {\displaystyle M} samples are used.
Note that solving for ฮป m {\displaystyle \lambda _{m}} will yield ambiguities, since only e ฮป m {\displaystyle e^{\lambda _{m}}} was solved for, and e ฮป m = e ฮป m + q 2 ฯ j {\displaystyle e^{\lambda _{m}}=e^{\lambda _{m}\,+\,q2\pi j}} for an integer q {\displaystyle q} . This leads to the same Nyquist sampling criteria that discrete Fourier transforms are subject to
| Im โก ( ฮป m ) | = | ฯ m | < ฯ ฮ t . {\displaystyle \left|\operatorname {Im} (\lambda _{m})\right|=\left|\omega _{m}\right|<{\frac {\pi }{\Delta _{t}}}.}
See also
โข Computation of Prony decomposition using Autoregression analysis
โข Application of Prony decomposition in Time-frequency analysis
Notes
References
โข citerefcarrieremoses1992Carriere, R.; Moses, R. L. (1992). "High resolution radar target modeling using a modified Prony estimator". IEEE Transactions on Antennas and Propagation. 40: 13โ18. doi:10.1109/8.123348.
โข citerefslyusar1998Slyusar, V. I. (1998). "Interpretation of the Proni method for solving long-range problems" (PDF). Radioelectronics and Communications Systems. 41 (1): 35โ39.